#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <pthread.h>

#define NTHREAD 100

double S, E, r, sigma, T, mean, *trials, *soma, *stddev, *t;
int M, tam;
pthread_t threads[NTHREAD];

double gauss_rand()
{
    static int    have = 0;
    static double x2;
    double fac, y1, y2, x1;

    if ( have == 1 )  /* already one available ? */
    {
        have = 0;
        return x2;
    }
    else
    {
        /* generate a pair of random variables */
        /* srand(time(NULL)); */
        y1 = (double)rand() / RAND_MAX ;
        y2 = (double)rand() / RAND_MAX ;
        fac = sqrt(-2*log(y1));
        have = 1;
        x1 = fac*sin(2*M_PI*y2); /* x1 and x2 are now gaussian */
        x2 = fac*cos(2*M_PI*y2); /* so store one */
        return x1;               /* and return the other. */
    }
}

void *calcula_desvio(void *parametro){

    int num = *((int *)parametro);

    int i;
    for(i=num*tam ;i< (num*tam) + tam; i++)
        stddev[num] += pow(trials[i] - mean,2);

    pthread_exit(0);
}


void *calcula_bs(void *parametro){

    int num = *((int *)parametro);
    int i;
    double max;

    for(i=num*tam ;i< (num*tam) + tam; i++){

        t[num] = (S*exp(((r-(0.5*sigma*sigma))*T) + (sigma*(pow(T,0.5)*(gauss_rand())))));

        //printf("%lf ",t);

        max = t[num] - E;

        if(max < 0)
            max = 0;

        trials[i] = ((exp((-1*r)*T))*max);
        soma[num] += trials[i];
        //printf("%lf ", trials[i]);
    }

    pthread_exit(0);
}

int blackScholesPar(void){

    int i;
    srand(time(NULL));

    scanf("%lf",&S);
    scanf("%lf",&E);
    scanf("%lf",&r);
    scanf("%lf",&sigma);
    scanf("%lf",&T);
    scanf("%d",&M);

    trials = (double*) malloc(sizeof(double)*M);
    soma = (double*) malloc(sizeof(double)*NTHREAD);
    t = (double*) malloc(sizeof(double)*NTHREAD);

    double soma_total = 0, stddev_total = 0;

    tam = M/NTHREAD;

    for (i = 0; i< NTHREAD ; i++){
                //Cria as threads do vetor com a funcao de calcular o Black Scholes
                pthread_create (&threads[i], NULL, calcula_bs, &i);
    }

    for (i = 0; i< NTHREAD ; i++){ //Espera todas as threads concluirem
                pthread_join (threads[i], NULL);
                soma_total += soma[i];
    }

    printf("SOMA TOTAL: %lf\n",soma_total);
    mean = soma_total/M;
    printf("MEAN: %lf\n",mean);

    stddev = (double*) malloc(sizeof(double)*M);

    for (i = 0; i< NTHREAD ; i++){
                //Cria as threads do vetor com a funcao de calcular o Black Scholes
                pthread_create (&threads[i], NULL, calcula_desvio, &i);
    }

    for (i = 0; i< NTHREAD ; i++){ //Espera todas as threads concluirem
                pthread_join (threads[i], NULL);
                stddev_total = stddev_total + stddev[i];
    }

    printf("STDDEV_TOTAL: %lf",stddev_total);

    stddev_total = pow((stddev_total/M-1),0.5);

    double confwidth = ((1.96*stddev_total)/(pow(M,0.5)));

    printf("confmin: %lf \n confmax: %lf \n", mean - confwidth, mean + confwidth);
}

